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Abstract 

We propose a method to reduce the relaxation time towards equilibrium in stochastic sampling 
of complex energy landscapes in statistical systems with discrete degrees of freedom by generaliz- 
ing the platform previously developed for continuous systems. The method starts from a master 
equation, in contrast to the Fokker-Planck equation for the continuous case. The master equation 
is transformed into an imaginary-time Schrodinger equation. The Hamiltonian of the Schrodinger 
equation is modified by adding a projector to its known ground state. We show how this transfor- 
mation decreases the relaxation time and propose a way to use it to accelerate simulated annealing 
for optimization problems. We implement our method in a simplified kinetic Monte Carlo scheme 
and show an acceleration by an order of magnitude in simulated annealing of the symmetric trav- 
eling salesman problem. Comparisons of simulated annealing are made with the exchange Monte 
Carlo algorithm for the three-dimensional Ising spin glass. Our implementation can be seen as a 
step toward accelerating the stochastic sampling of generic systems with complex landscapes and 
long equilibration times. 
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I. INTRODUCTION 



Numerical sampling of rugged energy landscapes is notoriously difficult [lj. Transition 
rates between two states are exponential functions of the energy barrier between them di- 
vided by the temperature. The sampling of complex systems is a daunting task because 
there are many states of comparable energies separated by large barriers. One of the most 
widespread sampling methods at finite temperatures is the Monte Carlo method, where one 
follows one or more trajectories of a virtual Brownian particle as it moves through the state 
space. To sample a system at a given temperature, there exists a plethora of approaches, and 
amongst many others particularly notable are the standard (Metropolis) [2] and the kinetic 
Monte Carlo [3] methods. At low temperatures these stochastic schemes tend to take a long 
time before a satisfactory result is reached. If one is interested in the behavior of a specific 
problem in the low-temperature limit, the common method is simulated annealing Jl]. If, 
however, one samples systems where below a certain temperature the state space splits into 
regions which are separated by huge barriers, like for example in spin glass systems, simu- 
lated annealing does not always lead to satisfactory results. Examples of accepted solutions 
to this problem include the exchange Monte Carlo method [5] and the population annealing 
[0] . Although good results are obtainable for most systems using these techniques, relaxation 
time still stays a crucial factor. For the protein folding problem [7] for example, the longest 
simulations available on present day computers are still far from the equilibrium distribu- 
tions. A possible step in resolving the problem of over-long relaxation times was proposed in 
Ref. [8] where a method to accelerate the sampling of continuous systems was introduced. 
The basic idea is to rewrite the Fokker-Planck equation, which describes the time evolution 
of the probability distribution for continuous systems, into an imaginary-time Schrodinger 
equation, for which one artificially introduces an energy gap between the ground state and 
the first excited state. Then the relaxation time, which is proportional to the inverse of the 
energy gap, is reduced. 

In the present paper, we extend the idea of Ref. [S] in Sec. [n] to discrete systems and 



implement it into a stochastic sampling scheme. Next, in Sec. Ill we analyze how our 
method can be used to improve the performance of simulated annealing of the traveling 
salesman problem [9] and the three-dimensional Ising spin glass. In the case of the traveling 
salesman problem we find that the simulated annealing is significantly accelerated and the 



modified sampling finds the approximately optimal solution much faster than unmodified 
simulations. Our method also leads to improvements for the three-dimensional Gaussian 
Ising spin glass. Sec. [TV] concludes this paper. 



II. ACCELERATED SAMPLING 

In this section we extend idea of Ref. [8] to discrete systems and discuss its implementa- 
tion into the kinetic Monte Carlo algorithm. We furthermore introduce a simplified version 
of the kinetic Monte Carlo scheme to speed up the calculations and save computational 
resources. 



A. Widening the Gap 



The master equation whose transition rates w a b fulfill the detailed balance condition 
reads, 

y>«^-yw« = e-w-'u) . (i) 



dPa 

dt 



b b 

If we set P a = f a Qa, with f a = e~^ Ea f 2 / \/~Z , where Z = J2b ex P(~fiEb) is the partition 
function of the system and E a is the energy of state a, we get 



dQa _ sr~^ fb 



^2 y W abQb - ^2 W baQa = ~ ^ H ab Q b . (2) 
b /a b b 

Since H a b, as defined in Eq. ([2]), is a real, symmetric matrix, we call it a Hamiltonian and 
thus, Eq. @ may be regarded as an imaginary-time Schrodinger equation. It has a zero 
eigenvalue with eigenvector fb, ^2 b H a t,fb = 0. This follows directly from the definitions of 
fb and H a b and the detailed balance condition in Eq. ([TJ. The lowest eigenvalue of H a b is 
therefore zero as guaranteed by the Perron-Frobenius theorem. 

Following the idea in Ref. [5j, we make the transformation H ab — > H ab + X(5 a b — -P° b ), 
where P® b is the projector to the state of zero eigenvalue, which is expected to shorten the 
relaxation time toward equilibrium. The spectrum of H ab is then shifted by A for all states 
except the ground state. It is easy to see that the matrix elements of the projector are 
Pa b = fafb, since Y.b P abfo = faJ2bfb = fa Efe exp(-/3E b )/Z = f a and the eigenvectors of 
H a b corresponding to different eigenvalues are orthogonal. Therefore, the imaginary-time 



dQg 

dt 



Schrodinger equation with the modified Hamiltonian is 

C (j™ ab + X^j Q b -Qa\X + ^2 Vbaj ■ (3) 

The solutions of this equation are, 

Qa(t) = Qi 0) + g«e^+ A > + Qi 2 )e"^ 2+A ) + • • • (4) 



where the {Q^} and {e n } are the eigenvectors and corresponding eigenvalues of H ab , re- 
spectively. 

Making use of the relation P a = f a Q a , we can translate Eq. ^ into the modified master 
equation 

( W <* + X< —^j Pb-PaU + Ys^j ( 5 ) 



dt 

b 



and deduce its solution from Eq. Q as 



p -l3E a 

Pa (t) = e —— + pW e -'te+A} + pf e -*(e 2 +A) + . . . (6) 

Zj 



For large times this decays to 



e -PE a 



Pa = (7) 

i.e. the Boltzmann weight, as expected. However, as can be seen from Eq. (|6]), the relaxation 
is faster than the case where A is absent. 



B. Implementation 

In Ref. [Bj, a similar idea was tested for continuous systems using a diffusion Monte 
Carlo calculation. In the present discrete straightforward implementation of our 

method is through the kinetic Monte Carlo algorithm j3]. Let us first describe this Monte 
Carlo method on the original master equation ([TJ) in order to make is clear what parts need 
modifications to accommodate the A-term in Eq. ([3]). The idea is to try to generate time 
'trajectories' of the system among its various available states in such a way as to satisfy 
the master equation. Assume the system is in a given state a. The rate (or probability per 
unit time) at which the system will escape from a to any available state b is given by the 
second term on the right-hand side of Eq. (jTJ). In other words, it is equal to T a = J2b w ba- 
Therefore, the probability distribution of the escape-time from a is given by the Poisson 
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distribution P^ sc (t) = T a e~ Tat . The time at which the system will leave state a can thus be 
drawn from this distribution. The probability for the system to go from a to a state b is 
obviously given by the ratio Wb a /T a . 

The practical implementation of the algorithm goes as follows. Assume that the system 
is in state a at time t n : 

• The system will make a random transition out of state a at a time t n+ \ = t n + r where 
t is drawn from the distribution P^ sc (r) = r a e _raT . In practice, one draws a number 
r uniformly distributed between and 1 and takes t = log(— r). 

• The state b to which the system will make the transition is chosen with probability 
Wb a /T a . A simple way to do this is to draw all individual probabilities consecutively 
until eventually they fill up the interval [0, 1] . Then one draws a uniform random 
number r\ between and 1, and the state b to which the system jumps is the one 
indexed by the transition probability which corresponds to T\ on the interval [0,1]. 

By generating many such trajectories, one generates probability distributions P a (t) which 
stochastically satisfy the master equation. Let us note that all trajectories generated this 
way are statistically independent, and can thus be used to compute averages. 

We now turn to an implementation of the kinetic Monte Carlo when one introduces the 
parameter A. Let us first expand unity as 



where the sum is over the states {c(a)} accessible from a given state a. If all states are 
accessible in principle, which is the case if the system has no intrinsic dynamics, like the 
Ising model, then we restrict {c(a)} to some subset depending on a, for example the nearest 
neighbors. Next, insert Eq. (M) into the outgoing part of the master equation, to obtain 



Note that, even if the physical considerations do not allow non-local moves, we can define 
transition probabilities w^ a = Xe~" Ec /Z a for such moves and thus incorporate non-local 
states into the list {c(ct)} in a very straightforward manner. 
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The introduction of the term Z a in Eq. (|8]) and Eq. Q makes the present implementation 
unsuitable to reproduce finite-temperature properties. To show this fact, let us first recall 
that a frequently-used transition probability is the heat bath (or Glauber) method 

which trivially fulfills the detailed balance condition. However, the modified master equation 
is not symmetric in its outgoing and incoming parts. The kinetic Monte Carlo algorithm 
uses only the outgoing part for sampling, and therefore, adding the A terms tilts the detailed 
balance in favor of the energetically lower lying states. To see this, we write the detailed 
balance condition as 



< m - + x Z„ 

w ac , x exp(-/3£ a ) 



(11) 



" ir„< I A 



Inserting the transition rate ( 10 ) into the above equation, we get 



wi = e-^(l + A[e-^ + e-^]/Z a ) = 0(Ec . Ea) Z c Z a + \[e~^ + e~^} 

e-^(l + A[e-^ + e -f> E °]/Z e ) Z a Z c + A[e~^» + e~P E °] ' 1 ' 

We see clearly that the detailed balance in its conventional form is not satisfied. If we 
recall that Z k = Yli m {k) e~^ Em , which represents the sum of Boltzmann factors of the subset 
{m(k)} of states accessible from state k, then we can conclude that if, Z c < Z a , there are 
more energetically lower lying states available from a than from c. This would suggest that 
in this implementation (Eq. (|9])) the addition of A indeed tilts the detailed balance in favor of 
states from where more lower lying states are accessible. This changes the finite temperature 
values of physical observables when compared to calculations made with A = at the same 
temperature. Nevertheless, when we are interested in the ground state solution, we do not 
have to worry about such finite temperature differences. 

There may be other implementations that do not use state- or temperature-dependent 
renormalizations of A and thus allow us to keep detailed balance. However, we reserve the 
finite temperature case for future studies. 



Using the language and example of the traveling salesman problem (see Sec. Ill A) we 
now introduce a simplified version of the kinetic Monte Carlo scheme. Sampling all nearest 
neighbors of a given tour, as needed for the kinetic Monte Carlo, is neither efficient nor 
feasible. At every Monte Carlo step we would have to calculate N(N — l)/2 transition rates, 



where N is the number of cities in the map. For example, if N 1000 we would have 
to make about 5 ■ 10 5 calculations. Therefore, we choose just a certain number of nearest 
neighbors and generate another number of non-local states to build our list {c(a)}. This is 
possible since the traveling salesman problem has no intrinsic dynamics, and thus, all other 
states are accessible from any given state. By taking a set of nearest neighbor states and 
a set of non-local states as possible jumps available, we sample in effect the local structure 
as well as 'far' away states and can in this fashion overcome large barriers. We can either 
choose to sample less local or non-local states than in the case of the kinetic Monte Carlo. 
In such a way the computational cost is reduced. 
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FIG. 1. (Color online) a) Simulated annealing of the traveling salesman map 'gr229' with 229 cities. 
Comparison of the full kinetic Monte Carlo (K) and the simplified (H) method, b) Investigation 
of the simplified method with v non-local and no local trials. See the text for details. 



Before performing full-scale computations, we run preliminary simulations to check if the 
simplified method actually works and investigate what values of various parameters are to be 
used in practical calculations. In Fig. [T] a) we compare the full kinetic Monte Carlo (labeled 
K) and the simplified method (labeled H). We perform simulated annealing calculations of a 
small traveling salesman map, 'gr229' [9J, with 229 cities. The annealing schedule is chosen 
to be step-wise growing (see Fig. [4] for comparison), starting from an already low j3$ = 40, 
the inverse temperature is increased by 5(3 every 2000 Monte Carlo steps. For the simplified 
method, we chose 230 local states and 230 non-local states per Monte Carlo step at random 
to build the list {c(a)}, and sample all nearest neighbors plus 230 non-local states for the 
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full kinetic Monte Carlo. The simplified method outperforms the full kinetic Monte Carlo, 
as the simplified method finds a shorter tour within the investigated time window. 

In Fig. [l]b) we investigate the parameter dependence of the simplified method. We look 
at the case where we allow no nearest neighbor hops and sample v non-local states per Monte 
Carlo step to see the effect of adding A to the transition rates while taking too few states to 
build the list {c(a)}. We see that for v = 25 and very small A = 0.0001 (red line) there is a 
slight visible deviation to shorter tours from the simulation with A = (black line), while a 
larger value of A = 1 (green squares) affects the sampling in a negative way. As an extreme 
case we take only v — 2 non-local states and A = 1 (blue circles), and see that the system 
does not relax. The reason for this failure at A = 1, v = 2 is that, with adding A to the 
transition rates, jumps are facilitated. Since we are looking at only two possible jumps, at 
some point both transitions will become approximately equally likely and the system fails 
to relax. We see that, if too few states to which jumps are possible are chosen, the sampling 
is influenced in an undesirable way. 

We use this simplified version of the kinetic Monte Carlo algorithm also for the Ising 
model, since for larger system sizes calculating the Hamiltonian and exponentiating it L 3 
times (there are L 3 nearest neighbors to a given state) becomes quickly very time-consuming. 
As a rule of thumb, we sample O(L), instead of L 3 , local states per Monte Carlo step. 

III. RESULTS 

In this section we begin with a short review of the models we used for our calculations to 
establish the terminology. Then we show the results of our simulations for different scenarios. 

A. The Models 

We applied our method to two models. First we treat instances of the symmetric traveling 
salesman problem. Given a set of coordinates M = {(x«, jji)} on a two-dimensional plane, the 
task is to find the shortest closed path, called a tour, connecting all points while traversing 
each point only once. Therefore, tours are ordered lists of the coordinates M, giving the rule 
in which order to visit the coordinates. We take the distance between two points % and j on a 
tour to be Euclidean d,E(i,j) = a/ (xj — Xi) 2 + (yj — yi) 2 so that the total length of the tour 
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is calculated as the sum of all segments E(r) = Y^^E(i,j)- To formulate this as a pseudo- 
physical problem, we identify the tours as the states, the tour length as the Hamiltonian and 
we choose nearest neighbor hoppings as the local dynamics. Nearest neighbors of a tour r 
are defined as tours s differing by the exchange of two points, i.e. having Hamming distance 
du of two to r, N(r) = {s|ci#(r, s) = 2}, see Fig. [2] 




FIG. 2. (Color online) Illustration of a nearest neighbor hop of the traveling salesman problem. 
The right tour differs from the left one by the exchange of the points indicated. 

With only nearest neighbor hoppings it is hardly possible to solve the traveling salesman 
problem to optimality [12J. We can, however, efficiently choose non-local hoppings [13], 
reversals or transport of tour segments. 

The second model we treat is the three-dimensional Ising model on the cubic lattice of 
linear size L, with random bonds. The energy of a given spin configuration S a is calculated 
by E a = — JijS^Sj — Sf , where the spins take values G { — 1, 1}. The bonds Jij 
are quenched random numbers drawn from the distribution P(Jij) = exp(— Jfj/2 J 2 ) / a/27tJ 2 , 
which we will call the Gaussian Ising spin glass, and h is an external field. We take only 
nearest neighbor jumps for the dynamics and use no non-local cluster flips [H], because 
we want to compare our method to the exchange Monte Carlo which gives satisfactory 
results with local sampling only. The nearest neighbors of a spin configuration S a are the 
configurations S b which differ from S a by a single spin flip, dn(S a , S&) = 1. 
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B. Simulated Annealing 

Simulated annealing [I] relies on one hand on a stochastic sampling scheme and on the 
other on an annealing schedule (3(n), the rule of increase of the inverse temperature with 
iterations. We first show on the wiggly harmonic potential discussed in Ref. [10] what 
the effect of introducing A is and then discuss the application of our method to simulated 
annealing. In this spirit we then apply simulated annealing to the traveling salesman problem 
and the three-dimensional Gaussian Ising spin glass. 

1. Wiggly Harmonic Potential 



Consider a one- dimensional system where a number of local minima is evenly distributed 
over the large basin U = x 2 /2, as in Fig. [3j The distance a between two neighboring minima 
is kept constant. The random walk of a point in this scenario is governed by the master 
equation describing the hopping between neighboring sites k — 1, k and k + 1: 



where B is the height of the barrier separating two minima and = a 2 ((k + l) 2 — k 2 )/2. 
When we take the continuum limit a <C 1, the coarse-grained master equation becomes the 
Fokker-Planck equation [TO] 




FIG. 3. The landscape of the wiggly harmonic potential. 




(13) 



dP(x, t) 
~~dt 



D 



d 2 P(x,t) 



+ /3D 



d[xP(x,t)] 



(14) 



dx 2 



dx 
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with diffusion constant D = e B/3 /a 2 . The fastest annealing schedule f3{t) that minimizes 
the average energy 

y(t) = / dxU{x)P{x,t) (15) 



is given by f3(t) = Int/B [10], which coincides with the generic bound for convergence to 
reach the global minimum [11J. With this schedule the average energy decays as y ~ 1/lnt. 
'Faster' schedules than this Int/B do not further minimize the average energy. 

We identify the cause for this inverse-log law for the wiggly harmonic potential as an 
instability in the associated imaginary-time Schrodinger equation. Then we will propose 
a way how a step-wise growing schedule, see Fig. |4j together with the considerations of 
Sec. |TT] can improve the performance. The average energy in this improved case scales as 
y(ti) = (Ei) ~ 1/5^, where i is the time step index and 5 is some constant, which implies 
y(t) ~ l/t. 



P 







time 



FIG. 4. Step-wise growing schedule. The inverse temperature is increased by a constant amount 
in regular intervals. 



To rewrite Eq. (14) as an imaginary-time Schrodinger equation, let us set P(x,t) 
e~^ u ^/ 2 7p(x,t) to get 



&ip(x, t) 
dt 



D d 2 ip(x,t) 



-D 



dx 2 

d 2 ijj(x, t) 
dx 2 



D^U"(x) 



(3U{x) 



V>0M) 



D 



(3 (3x 2 



4 



ip(x,t), 



(16) 



where $ — d/3/dt is the time-derivative of the schedule. We can draw on an analogy with 

/ .\V 2 

the quantum harmonic oscillator (h = 1) to have m = 1/2D, tu = D 1 ^ 2 \D(3 2 — $) ■ For 



a meaningful analysis we require that the separation of energy levels, which is proportional 
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to the frequency u, be greater than or at least equal to zero, 

D/3 2 - $ > 



(17) 



or 



3 /3 2 > a 2 /3. 



(18) 



This relation is asymptotically only fulfilled for functions which grow like In t or slower than 



that. To see this we first integrate Eq. (18) and have 



t > 



d/3 



a 2 e B ?. 



(19) 



We take the large (3 limit and ignore (a/(3) 2 compared to e BlS to get the usual Int/B > (3 
schedule restriction. For e~ B ^(f3 /a) 2 = (3 the energy levels coalesce and Eq. (16) becomes a 



free diffusion equation with time dependent, monotonically decreasing diffusion coefficient. 



Schedules which do not fulfill the relation (18) have an imaginary u and do not lead to a 
decaying solution. 

To circumvent this problem we employ the method proposed in [8] (see also Sec. [n]) for 
the accelerated sampling of Boltzmann distributions. If we choose the schedule like in Fig. 
El (3{t) = (3iQ(t — £j), where Q(t) is the Heaviside step function, then for intermediate 



U < t < t i+ i times, the Schrodinger equation (|16|) reads 
dtp{x,t) 



dt 



d 2 jj(x,t) 

OX 



Pi 



-X 



2 



ip{x,t) = Hii/j, 



(20) 



where (3i = Yl]=ofij an d = ^(A) = e B,3i /o 2 . Again, in correspondence with a harmonic 



oscillator, we obtain 



1 



rn 



2A 

uji = (3iDi 



E 



(») 



nu>i = nfiiDi = n(3id 



(21) 
(22) 
(23) 



and denote the eigenfunction of Hi in Eq. (|20|) to the eigenvalue as 



We now 



employ the sudden approximation [16], which uses the fact that, if the system changes too 
quickly, the wave function cannot follow and we can use the new Hamiltonian in the time 
evolution operator with the previous wave function as the initial condition. Let us assume 
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that a particle is in the ground state of the Hamiltonian at time step i — <Pl_.-i- Then just 
after the jump to the next time step i, the wave function of the particle can be described by 

j,(t) = e-^+'VS- (24) 
We insert unity as 1 = 1 0^X0^1 into the above equation to find 

4 = Y,e-^ iU «H$Wi)4 n) . (25) 



n 



Note that for n odd, the overlap vanishes since <f)f_i is an even function. To keep the 
probabilistic interpretation of e~^ u ^ 2 -ip, the normalization is chosen as 

exp (-— J iftdx = 1. (26) 



4 

Therefore, the normalized ground state solution of the Hamiltonian ifj_i reads 

m=(¥T°* (27) 

and the new wave function up to the slowest decaying term is 

^ = + m^e-*^^. (28) 

2 

For our approach to be sensible we have to wait again long enough for the decay of the 



excited state <p\ . However, the decay constant Tj = , as we learn from Eq. (23), is 

exponentially increasing with {3, and therefore, we will have to wait longer and longer as 
time proceeds for the system to decay. 

Now, we set Hi — > H \ + A(l — P,^), with P^ \ the projector to the ground state. Then 



Eq. (28) becomes 



1> = + m^e-^^<^. (29) 



(2) 

For any monotonically growing j3, E\ vanishes exponentially, so that the slowest decaying 



terms of Eq. (29 ) decay approximately as e A (*»+'). Therefore, the decay constant is bounded 



from below by Tj = 1/A and is thus independent of the index i. 
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Let us investigate what the above considerations mean for the average energy at a time 
t > ti, shortly before the next jump of (3(t) at £ i+1 , 

. x * — ( P^ 2 \ \/AA-i ,(o) 



Pi + Pi-i 
2 



VM^~i f,x 2 ( f3 t x 2 \ ( pA 1/2 ( (3x< 
' ax— exp — I I — exp 1 



A + A-i J 2 4 7 V2tt 

2 

ft \ 1/2 v 7 ^ 



27T 7 Pi + Pi-l ^3/2 

2 



A + A-i Pi 



(30) 



where we have used Eq. (27). For the schedule Pi = P + Si the prefactor becomes 
ygg^ _ 2^(gTgj(g+fcg^j)) ^ 2fe^g£j) 

A + ft-i ~ 2/9b + 2iW-ife ~ fc(2i-l) torz»l. (dlj 

2 

Thus the average energy scales as (E) ~ in contrast to (E 1 ) ~ B/ln(i), as is the case 
for the logarithmic schedule. 



2. Traveling Salesman Problem 

We now turn our attention to the traveling salesman problem. We use the 'gr666' and 
'ul060' data sets from the TSP-database [9] with 666 and 1060 cities, respectively. The 
inverse temperature is chosen to be stepwise growing with Monte Carlo steps. Starting 
from a base value Po = 40, the inverse temperature is increased by 5 = 3.5 after 3000 Monte 
Carlo steps, similarly to the plot in Fig. |4| For the simulation we used the simplified method 



described in Sec. |IIB| At each Monte Carlo step, we chose iV nearest neighbors and iV non- 
local states [13] randomly, where iV is the number of cities. The results of the simulations 
are shown in Fig. [5] a) for 'gr666' and b) 'ul060', where the tour length is plotted versus 
the logarithm of the Monte Carlo steps. The effect and advantage of using our method are 
clearly visible. Sampling for the same amount of iterations, we find much better solutions 
by using the transition rates modified by A, as defined in Eq. ([9]), than when using A = 0. 
Stated otherwise, we can achieve an acceleration by an order of magnitude to reach a given 
tour length. 
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a) b) 




Monte Carlo steps Monte Carlo steps 



FIG. 5. (Color online) Results of simulated annealing for the traveling salesman problem: a) 'gr666' 
and b) 'ul060' data sets with 666 and 1060 cities, respectively. 

3. Three-dimensional Gaussian Ising Spin Glass 




FIG. 6. (Color online) Comparison of the results obtained by simulated annealing (A) and exchange 
Monte Carlo (X) of the three-dimensional Ising spin glass, a) Magnetization per spin b) Edwards- 
Anderson order parameter. Note that the finite-temperature values of the present method (curves 
marked (A), A > 0) do not correspond to equilibrium properties. 



Next we treat the three-dimensional Gaussian Ising spin glass (see Sec. Ill A). It is clear 



that a naive simulated annealing, with local jumps only, fails to find a good low temperature 
solution of this model in a reasonable time. There are many energetically close states which 
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are separated by huge barriers, and thus the time needed to escape from a local minimum 
increases very quickly as the temperature decreases. Nevertheless, we would like to see 
how our method affects the annealing procedure. The temperature is lowered in 480 steps 
from T max = 5/J to T m ; n = 0.5/ J. The coupling constants of the Hamiltonian are drawn 
from a Gaussian distribution, with zero mean and variance J. The external field strength 
is chosen as h/J = 0.1. In Fig. [6] we compare the a) magnetization and b) Edwards- 
Anderson order parameter, q — /(-^ J2i ^i^^i^) 2 ^ °f simulated annealing (A) with results 
from exchange Monte Carlo calculations (X). Notice that the finite-temperature values of 
the present method, marked (A), do not represent the equilibrium properties for the reason 



discussed in Sec. |II B| The exchange Monte Carlo is usually expected to give good results 
although we have not checked equilibration conditions since our goal is to compare the 
performance of the methods under the same conditions on computational cost. 

The first observation from the data is that the introduction of the A term significantly 
improves the performance of simulated annealing at the lowest temperature. The values of 
physical quantities, M/N and q, have come close to those of the exchange Monte Carlo, the 
latter being a benchmark. Another notable fact is that the A term, at least for a small value, 
induces no perceptible change in the exchange Monte Carlo. Lastly, the very large value of 
A = 1000 yields close results to those of the exchange Monte Carlo at the lowest temperature. 
The results for A = 1000 are, nevertheless, still slightly away from the exchange Monte Carlo 
values. These facts clarify the usefulness as well as limits of the present method for this 
problem of Ising spin glass. 



IV. CONCLUSION 

Based on the idea of Ref. [8], we introduced a method to accelerate stochastic sampling of 
discrete, classical problems. Our method suggests a way to overcome the limits of standard 
simulated annealing. We tested our algorithm on the traveling salesman problem, where, 
in the framework used in this paper, we find the shortest tour an order of magnitude faster 
by using our method than in the conventional case. Simulated annealing of the three- 
dimensional Gaussian Ising spin glass is also accelerated. In this latter case, the performance 
of our method is relatively close to that of the exchange Monte Carlo. 

Throughout our investigation we used a simplified version of the full kinetic Monte Carlo 
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algorithm to reduce the computational cost of the sampling at each Monte Carlo step. 
This simplified method outperforms the full kinetic Monte Carlo when we are faced with 
a plethora of accessible states, but when the choices are limited, adding A tends to have 
undesirable effects if it is not chosen accordingly. However, in its present form the algorithm 
is useful only in the search for very low temperature solutions. 

In conclusion, the present method would be a useful alternative of simple simulated an- 
nealing for optimization problems. The relatively straightforward implementation using ki- 
netic Monte Carlo and non-local moves would make it a method of choice for some purposes, 
especially where computational cost is a factor. 
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